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Estimation of Ambiguity Functions With Limited 

Spread 

Heidi Hindberg & Sofia C. Olhede 
Abstract 

This paper proposes a new estimation procedure for the ambiguity function of a non- stationary time series. 
The stochastic properties of the empirical ambiguity function calculated from a single sample in time are derived. 
Different thresholding procedures are introduced for the estimation of the ambiguity function. Such estimation 
methods are suitable if the ambiguity function is only non-negligible in a limited region of the ambiguity plane. 
The thresholds of the procedures are formally derived for each point in the plane, and methods for the estimation 
of nuisance parameters that the thresholds depend on are proposed. The estimation method is tested on several 
signals, and reductions in mean square error when estimating the ambiguity function by factors of over a hundred 
are obtained. An estimator of the spread of the ambiguity function is proposed. 

Index Terms 

Ambiguity function, estimation, harmonizable process, non- stationary process, thresholding, underspread process. 



I. INTRODUCTION 

We propose a new estimation procedure for the Ambiguity Function (AF) of a non- stationary Gaussian process. The suitable 
definition of a 'locally stationary process', and the development of inference methods for a given sample from a non- stationary 
process is still an open question, see | IJ. The AF forms an essential component to this task, as it provides a characterisation of 
dependency between a given time series and its translates in time and frequency. The AF of a non- stationary process can often 
be assumed to be mainly limited in support to a small region of the ambiguity plane, and this corresponds to correlation in the 
process of interest being limited in support. If the support is additionally centred at time and frequency lags of zero, then the 
process is underspread, see Matz and Hlawatsch O. An important class of non- stationary processes, namely semi- stationary 
processes |3|, exhibit limited essential spread in the ambiguity plane. Non- stationary processes can also be constructed from 
the viewpoint of time- variant linear filtering. Pfander and Walnut [4] have shown that if the spreading of a filtering operator 
is sufficiently limited then the operator may be identified from its measurements. It is not necessary in this case to assume 
that the spread is centered at the origin |4 |. The key to the estimation or characterization of the generating mechanism of 
non- stationary processes, in all of these cases, is the compression of the AF of the process. 

The AF has also been popularly used in radar and sonar applications 0, where the echo of a known signal is recorded. 
The delay and frequency shift of the echo can be determined from the AF of the signal |6|. In such applications the emitted 
signal is deterministic, even if the echo has been contaminated by noise. By using the compression of the emitted signal in 
the AF domain, the properties of the process of reflection can be determined, see Ma and Goh (71. The estimation of the AF 
of a non-zero mean signal immersed in noise is of interest. Unlike Ma and Goh, we shall not assume the compression of the 
AF is known, but propose an estimator of the AF, suitable if the AF satisfies some (unknown) compression constraints. 

Given the importance of the AF it is surprising to find that existing methods for its estimation are quite naive. Methods are 
usually based on calculating the EMpirical AF (EMAF), or some smoothed version of the EMAF. It is clearly unpalatable to 
implement a smoothing procedure on the EMAF, as the compression of the EMAF should be preserved, or even preferably 
increased, by any proposed estimation procedure. In this spirit Jachan et al. | 8 1 have used shrinkage methods to estimate the 
EMAF with an assumed knowledge of its support, using a multiplier that attenuates larger local time and frequency lags. Often 
the support of the AF is not known, and arbitrary shrinkage at larger lags will not be an admissible estimation procedure. 

Given the compression or sparsity of the AF, shrinkage or threshold estimators are a natural choice. In this paper we propose 
a set of threshold estimators for the AF. Thresholding has been implemented before in the global time-frequency plane for 
wavelet spectra, or evolutionary spectra, see Fryzlewicz and Nason |9|, or von Sachs and Schneider I.1QJ . The estimation 
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problem in this case is quite different to estimating the AF, as we cannot assume the raw wavelet/evolutionary spectra exhibit 
compression. Hedges and Suter |11| calculated the spread of the AF, based on exceedances of the magnitude of the EMAF 
over a (fixed and user tuned) threshold in a given direction in the local time-frequency plane. The stochastic properties of the 
EMAF were not treated by Hedges and Suter, and instead their investigation focused on some effects from different boundary 
treatments. 

To be able to select a suitable threshold procedure at each local time-frequency point we must establish the distribution of 
the EMAF, and this requires modelling the observed process. We discuss the classes of non- stationary processes that we intend 
to treat in Section |II-A| and initial method of moments estimators of second order structure in Section |II-B| In Section [In] we 
determine the first and second order properties of the EMAF. For strictly underspread processes satisfying constraints on their 
degree of uniformity we show that the EMAF is asymptotically Gaussian proper, if calculated for local frequency and time 
lags outside the support of the AF. 



We introduce the proposed estimation procedure in Section IV for deterministic signals immersed in noise and stochastic 



processes. In the first estimation procedure we threshold the EMAF based on a variance estimated from the entire plane, 
this yielding the Thresholded EMAF (TEAF). For stochastic processes, we also introduce a Local TEAF (LTEAF), where 
we estimate the variance locally in regions of the ambiguity plane. Such procedures are suitable as smooth time-frequency 
spectra will correspond to decaying ambiguity functions and are a natural extension to treating wavelet coefficients differently 
depending on the level of the wavelet transform, so called level-dependent thresholding, see e.g. Johnstone and Silverman 1 12|. 
The EMAF of deterministic signals immersed in analytic white noise will be biased. We propose a method to remove this bias 
prior to local thresholding, thus obtaining the Bias-corrected LTEAF (B LTEAF). We define an estimator of the total spread of 



a signal in the ambiguity plane in Section |IV-A[ to numerically measure the support of the signal in the ambiguity plane. The 
proposed TEAFs will produce sparse approximations to processes that are not strictly underspread. Whenever the EMAF is 
small in comparison to the variance of the EMAF, the AF will be estimated as zero. The AF can provide essential information 
as to the natural bandwidth of a process via the spread of the TEAF. 

In Section [V| we provide simulation studies of the proposed estimators. The Mean Square Errors (MSE) of our procedures 
reduce in many cases to around a hundredth of the MSE of the EMAF. Plots of single realisations show the accuracy of the 
threshold estimators, based on a single sample. The proposed methodology thus introduces a new method of characterising the 
features of structured non- stationary processes with high accuracy. 

II. Second Order Structure Descriptions 

A. Modelling the Observations 

We assume that the process under consideration, X[t], is an analytic, zero-mean, discrete-time, Gaussian harmonizable 
random process. Note that [•] is used to indicate a discrete argument, and (•) is used for a continuous argument. As X[t] is 
harmonizable it admits the spectral representation 1 13] of 

X[t\= / e^^^'UX{f), tez, (1) 
Jo 

where is a Gaussian increment process, and X[t] G C. We note that has a correlation structure specified by 

the dual-frequency spectrum Sxx* /) du df = E + The AF of X[t] is defined as a Fourier Transform 

(FT) pair with the dual-frequency spectrum in variable /. The AF forms an FT pair with the dual-time second moment of the 
process, or Mxx* r] = E {X[t]X*[t — r]}, for t and r G Z, now in variable t. We refer to 1 141 . [151 for a more detailed 
exposition of the AF and forms of AFs of common processes. 

Estimation of the AF of a real-valued process is known to display artifacts from aliasing and interference from negative 
frequencies 1 16|. For this reason we consider analytic processes solely. Note that since we are working with the analytic signal 
of a non- stationary process, it may (in general) be improper ifTTl and thus may exhibit complementary correlation |18|. We 
leave the estimation of the complimentary AF outside the scope of this paper, noting that our proposed methods can in a 
straight forward fashion be extended to include the estimation of such objects. Furthermore, we will in this paper assume that 
we are working with proper processes, and any terms that will only be non-zero for improper processes will be omitted. 



B. Initial Estimation 



We observe a finite sample of a realization of X[t]. The Sample AF (SAF) of a sampled deterministic signal, {g[t]}^^Q^ , 
is given by 



Ar-l+min(0,r) 
t=max(0,r) 



1 1 

2' 2 



_(7V-l),...,(iV-l)]. 



(2) 
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The properties of Agg*{iy^r] follow directly from properties of AFs of deterministic sequences, see |[T9l , 1201 . The effects 
of discretization on deterministic structure, i.e. trying to infer properties of g{t) from g[t] = g{tAt) for some appropriate 
sampling period At > 0, will also be left outside the scope of this article [211 . 1. 22 1 . The Empirical Second Moment (ESM) 
of a stochastic process based on the realization x[t],t = 0, . . . , — 1 is given by Mxx* t] = x[t]x*[t — r], which has an 
expectation E |Mxx* [t, t] | = Mxx- [t, r] and a variance var |Mxx* [t, t] | = Mxx- [t, 0]M^^. [t - r, 0]. Note that the 
ESM, its expected value and variance is defined as above for < t < — 1 and t — {N — 1) < r < and they are zero 
otherwise. Clearly the variance of Mxx* ^] is very large and the estimator should not be put to direct use. We form the 
EMAF of X[t] by 

Ar-l+min(0,r) Ar-l+min(0,r) 

Axx'{v,t]= Mxx'[t,r]e-^''"'' = E x[t]x*[t - rje-^^^'^K (3) 

^=max(0,r) t=max(0,r) 

The EMAF corresponds to an estimator of the AF of X[t]. However, it is not a suitable estimator of Axx^i^j^], as 
its variance will be extremely large. When implementing ([3| we have chosen zero-padding as our choice of extending the 
estimated covariance matrix. We found that periodic (circular) extension created unpalatable mixing effects of a fundamentally 
non- stationary object. Edge-effects are discussed more carefully in 1.11 J . The best choice of edge treatment will depend on the 
statistical and deterministic properties of the chosen EMAF. 

III. Properties of the Empirical Ambiguity Function 

A. Representation and Moments of the EMAF 

We shall construct estimators of the AF, based on the EMAF. The EMAF of an analytic zero-mean harmonizable process 
X[t] admits the representation of 

Jo Jo 

where D^if) — sin(7r/A/')/ sin(7r/). We respectively denote /i^l^, t], cf\{-u^t] and rA(^, t] as the mean, variance and 
relation of the EMAF. We denote an analytic process corresponding to a real-valued stationary white process as an analytic 
white process, but keep in mind that this process is not actually white. Thus, an analytic white process is a stationary process 
that has a power spectral density that is constant for nonnegative frequencies, and zero for negative frequencies. 

Proposition 3.1: Moments of the EMAF for noisy Deterministic Signals 
For a process that is the sum of a deterministic analytic signal and analytic white noise with variance a^, X[t] = g[t]-\- W[t], 
the expected value, variance and relation of the EMAF take the form 

fiAit^, t] = Agg, {v, r] + ^e-^---(^+--i)D^,_|,| {v)e^^^'\mc{T 12) (4) 

^^(^^, r] = al^h{v, t] + a^wH-^, -t] + cj^w^N - |r|)(l/2 - \v\) + O (1) (5) 
rA{v, t] = -<J^\v\{N - \t\)L (N - T, v) e-^^^''^''-^'^smc{2\v\T) + 2al.h'{v, t]+0 (1) 

/oo />l/2 — max(0,z/) 

sinc(/)sinc(/ + 2(iV-|T|Md/, h{v,T]= |G(/,T]|'d/ 
-CO ^max( — i/,0) 

/>l/2+min(0,-i/) 

h'{v,T] = / G*(/,r]G(/ + -r]e^^^^^-''^<^>^^df 

J max(0, — z^) 

Ar-l-|r| 

where sinc(x) = sm{iTx) / {irx) and G(/, r] = XI ^[^ ~ kl^C'^ < 0)]e~-^^^-^*. The indicator function is defined as 

/(r < 0) = 1 if r < and /(r < 0) = otherwise. 

Proof: See appendix [A[ ■ 
Proposition 3.2: Moments of the EMAF for a Stationary Process 
For a zero-mean, analytic stationary stochastic process X[t] the expected value, variance and relation of the EMAF take the 
form 

^iA{v,T] = D^,_|,|(i.)e-^-'^(^+--i)Mxx.M (6) 
a\{v,T] = {N-\T\){l/2-\p\)Axx'{-'y,0]+O{l) (7) 
rA{^,T] = e-^^^''^^+^-'\N-\T\){l/2-\,y\)L{N-\T\,v)Axx'{'y,T]+0{l) (8) 
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where Mxx* [^] is the autocorrelation sequence of the process, Sxx* (/) is the spectral density of the process, and 

/>l/2+min(zy,0) _ _ 

Axx- iy, r] = / Sxx' if - i^)Sxx' Uy^^^^ df/ {1/2 - \v\). (9) 

J max(0, — zv) 

Proof: See appendix |B] ■ 
The AF of a stationary process will be nonzero only on the stationary manifold = 0, and takes the form Axx* ^] = 
Mxx* For a fixed value of A^, I)7v-|r|(^) does not correspond to a delta function in v. Ideally however the estimated 

EMAF should not exhibit spreading in u. 

Proposition 3.3: Moments of the EMAF for a Uniformly Modulated White Noise Process 
For a zero-mean, analytic process X[t] corresponding to a real- valued uniformly modulated white noise process, the expected 
value, variance and relation of the EMAF take the form 

fiA{i^,r] = (l/2-|^|)Sxx*Me^'"(^/'-l^l)^sinc((l/2-|^|)r) + 0(l) (10) 

/l/2-\u\ ly^^ /fA|2 
I XX / I ^j2.fr^f^o{l) (11) 
-i/2+\u\ ^ - rl 



/»l/2+min(0,z^) /»l/2+min(0,i/) 
^max(0,i^) max(0,i^) 



where is the DFT of the time- varying variance of X[t], crj^lt]. 

Proof: See appendix [C] 

We note that the result of the integral in the expression for the variance decreases with and also with |r|. 



B. Distribution of the EMAF of an Underspread Process 

To determine suitable estimation procedures for the AF we need to derive the distribution of the EMAF. 

Theorem 3.1: Distribution of the EMAF of an underspread process 
Assume X[t] is the analytic process corresponding to a harmonizable real- valued zero-mean Gaussian process whose EMAF 
is strictly underspread. A strictly underspread process has an ambiguity function that is only non-zero for (z^, r] G P, where 
there exists some finite non-negative T and Vt such that V (Z [—Vt.,Vt\ x [— T, T]. Then the EMAF of X[t\ evaluated at [y^r] 
has a mean, variance and relation of 

1/2 1/2-/ 

/iA(z^,r] = ^xx*(z^^/)e^■'-^^e^-(^+--^)(-'--)i^^_,.|(z/'-z.)^z.'# (12) 



-/ 

A/'-l+min(0,r) T-1 



^=max(0,r) r' = -(T-l) 





'TV' 


^log 






T 



Then, with — > denoting convergence in law 1 23 1 



(13) 



while if < 1] and \r\ < T then the relation in noted in ( |A-8| ). Fix (i^, r) ^ V. Let r > and take 

Qiv-|r| [t] = e-^'2-^(^+^) {X[t + r]X* [t] - Mxx* + r, r]) , t = 0, . . . , N - 1 - r. (15) 
Assume that the triangular array {Qat-ItI t = 0, . . . , TV — 1 — r} is strongly mixing and satisfies for any e > 

^ N-T-l 

^- var{QAr-|r|M}^(|QAr-|r|MI >eaA(^,r]) ^0, (16) 

as A/" — |r| ^ oo, as well as the condition: 

^ N-T-l 

sup^- V Mxx*[t + r,0]Mxx*[t,0] < oo. (17) 



AT^ (0,1,0), (z/,r)^P, t>0, (18) 
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where J\f^ (0, 1, 0) denotes a complex Gaussian distribution with mean 0, variance 1 and relation 0. The result follows mutatis 
mutandis for r < 0. 

Proof: See appendix [P] ■ 
It is tempting to attempt to deduce the distributional results directly from the first and second order structure established in 
the first part of the theorem. In fact, if X[t] corresponds to a stationary process, where T and r are both even, the result 
follows directly. For more general classes of processes the result is slightly more involved, and a condition on the variance 
needs to be combined with a suitable mixing assumption, as above. Note that the CLT theorem that we use does not provide 
rates of convergence. Also, in some degenerate cases the joint distribution of \ Axx* t] > may not be asymptotically 
multivariate Gaussian. For ease of exposition, such cases are not treated here. Finally, whilst (asymptotically) retrieving the 



Gaussian distribution for the EMAF by Theorem 3.1 we fail to obtain simple interpretable forms for (j\{i'^t] and rA(i^, t]. 
For this reason, it is convenient to derive the first and second order structure directly from the modelling assumptions of some 
commonly used processes, as done in Section lllll 



IV. Estimation Procedure 

To determine the approximate distribution of the EMAF of a deterministic signal immersed in analytic white noise, we need 
to note the distributions of the four components of \A-\\ . Here, Agg* {u, r] is constant, whilst Aww* (^^ ^] is asymptotically 
Gaussian, and its form is determined by Theorem |3.1[ We note that since W[t] is Gaussian analytic white noise, and so 
from ( |A-3| ) and ( |A-6| ) 

Awg* r] + Agw* r] = J\f^ (O, {h{u, r] + h{-y, -r]) , 2(jlrh\u, r]) . 



From (A-2) and (A-5) we can note that Awg*{i^jT] + Agw*{i^^T] is independent of Aww^iy^^]- Combining (A-1) and 
Theorem 3.1 as well as using proposition 3.1 it follows for {f^r) ^ (0,0) 



(0,l,r^(^,r]M(^,r]) 



(19) 
(20) 



Hence we may note \iia{i^-,t]\^ <^ cf\{u^t] in most of the ambiguity plane. Furthermore it follows that 



(j\{v, t] 



VAiy. r]| 
(7\{y, r] 



Ul- 



0(1), 



(21) 



where IJ\ and IJ2 are independent standard Gaussian random variables. Thus for such points in the ambiguity plane it follows 

^ 2 

that Axx* t] — is a weighted sum of Xi's. These are equally weighted if rA{^, r] = 0. For most points in the 

ambiguity plane this statement will be roughly valid as apparent from proposition 3.1 Then if (z^, r] ^ V, where V denotes 
the support region of the ambiguity function of the real part of X[t], we may note that 



Axx*{j^,r] - iiAij^.r] 



;X2 



O 



N-\t\ 



-0(1), 



(22) 



A suitable estimation procedure for such random variables with cr\(v^T\ known a priori, is then the following thresholding 
procedure, for some given threshold A^, 




A-xx' {v, t] 



<\^a\{v,T] 



(23) 



If we normalize 



via dividing the sequence, element by element, by cr^(z/, r], then we retrieve a set of 

correlated positive random variables. For any such collection, we may note from 1241 . that using a threshold A^^(C) = 
2\og{NxWog{Nx)]^) for C > 1 is suitable. In our example Nx is chosen to be twice the number of observations, as we 
threshold the ambiguity function frequency by frequency, across the total collection of all time lags. Olhede |24|[p. 1529] has 
calculated the risk of this non-linear estimator for sums of unequally weighted Xi's. 
Of course (j\{v^t] is not known. We standardize the EMAF by 



A 



XX 



^{N-\T\){l/2-\v\r 



(24) 
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and note that 



var 



^4 



2 H^^^ 



h{- 



-r] 



(7V-|r|)(l/2-|H) 



O 



N ■ 



(25) 



We have assumed that g[t] G 
at high values of and |r| 



£2, and so \h{i'^r]\ = 0(1). Thus a 



2 h{u,T]-\-h{-u,-T] 



= 0(1), where the form of h{iy,T] 



W iN-\T\)il/2-\u\) 

ensure that the statement is still true for such values. We shall use the Median Absolute 



Deviation (MAD) estimator of the variance. MAD has been used for estimating the scale of correlated data before, see 
|[T2ll . The MAD estimator will need an adjustment factor that is different for ^xl random variables compared to Xi random 
variables. We note the median of a ^X2 In [2]. We define an estimator of the analytic white noise variance for any region 
{(v, t]}=Mc [-1/2, 1/2] X {-{N - 1), . . . , (iV - 1)} by 



median 



iM) = 



ln[2] 



(26) 



The imprecision of this procedure will depend on the lack of compression of the representation of in the ambiguity 

domain. We note that MAD has a breakdown point of 50 %, and so with quite severe contamination the estimator will still 
be useful, if somewhat inefficient. We then take 



al{u,r\M]=dUM){N-\r\){l/2-\u\). 



(27) 



A suitable threshold procedure (now with an unknown (j\{i>^t]) simply corresponds to using (23) with the variance replaced 
by its estimated value from ( [27| ). If the variance is estimated from the entire plane, i.e., M = {—N/{2N)^ . . . , A^/(2A^)} x 
{-{N -1),...,{N - 1)}, we denote 1^^* {u, r] the Thresholded EMAF (TEAF). 



The EMAF of a deterministic signal immersed in analytic white noise is biased, see proposition 3.1 Given an estimator of 



cr^, we can remove this bias prior to thresholding. We define the bias-corrected EMAF by 



^XX 



'w 



{Mi)e-^'"''-^+^-^^DN-\r\ (i^) e-^'^^^/^sinc (t/2) . 



(28) 



The noise variance is estimated from some given region A4i, which is chosen to only include the rims of the ambiguity 
plane. We normalize A^x* ^] ^24). For smaller sample sizes, it may be unreasonable to assume that the contributions 
of h{iy^r] are negligible compared to \/ N — |r|. We define ^x(^, t] = cf\{u^t]/[{N — |r|)(l/2 — \i>\)\. This is explicitly 
subscripted by X, as it is only non-constant from contributions due to h(y^T\. We propose to segment the plane into regions, 
and implement the proposed procedure after taking a separate estimate of o\{y^T\ in each region. The optimal choice of 
region is given by (25). Since h(y^T\ is unknown, we propose to use a centre square around (0,0], and square annuli, see 
Fig. [ija), to estimate the local variance. This is based on assuming G\{y^T\ smooth and decaying from (0,0). In general we 
separate the ambiguity plane into K regions {A1/c}, so that '^\{y-, r] a^(A^/c], a constant for all (y^T\ e Mk- We estimate 
the variance in each region as 



median 



^xx 



ln[2] 



(29) 



We threshold A^^]^^{i'^r] according to (23), but now using cr^xi-^k) in each defined region instead of a^{A4) in (|27|, this 



yielding A^^*^(z/, r], or the Local Bias-corrected TEAF (LBTEAF). 

Next, we derive results for estimating the EMAF of an underspread zero-mean stochastic process. This is strictly speaking 
a correct treatment for processes satisfying the constraints of Theorem |3.1[ but we expect that the distributional result is valid 
under less constrained scenarios. The conditions are sufficient, but by no means necessary, for the distributional result to hold. 



We note from Theorem 3.1 that for most points in the ambiguity plane 

2 



A 



XX 



(i/,r] - /iA(^,r] 



(i 1 2 
= ^X2 



o 



N-\t\ 



0(1). 



(30) 



The question then arises of yet again estimating a\{v,T] appropriately. We note from Theorem 3.1 that a\{v,T) takes 
different forms depending on the spreading of the process. The more X[t] is like an analytic white noise process, the less 
variation will cr\{v, t] exhibit from the form of the variance of the EMAF of an analytic white noise process. We again define 



,(^,r]=<Ti(^,r]/[(iV-|r|)(l/2-|H)].Then 



(7V-|r|)(l/2-|H) 



'X 



O 



N-\t\ 



0(1) 



(31) 
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Fig. 1. (a) Partitioning of the ambiguity plane, (b) Radial lines over which the integrand in ([32| exhibits stationary phase. 



Using inverse FTs we note from Theorem 3.1 that the variance of the EMAF of a strictly underspread process can in fact be 
rewritten for r > 



T-l 



A*xx'i'^",T'y''^^-^-^^<-'''-''"^DM-ri'^' - v") dv' dv" + O ( lo; 

T-l M 



(32) 



E / e-^'^^''^'-'''-^\Axx'{i^',r']\Ui.' + o(\og 



0(1). 



Deriving this requires that Axx* t] is sufficiently smooth in v. If this is not the case, use (13). The variance of the EMAF 
corresponds to aggregating the total magnitude squared of the AF over the plane and implementing phase shifts. If the function 
is smooth in v with a stationary phase approximation to the integral we mainly pick up contributions on the line v' = vr' /r, 
which we later sum over r^ This is exemplified for some choices of u and r in Figure [TJb), where we plot the lines summed 
over. Despite this geometrical intuition, the form in (32) is not extremely informative. Equation (32) does indicate that cr^(z/, r] 



should be a smooth function of u and r, as we only integrate over V. We rely on the forms of propositions |3.2| and |3.3| to 
give more precise understanding of the variance of the EMAF in these special cases. We note from these set of results that the 
variance is smooth in v and r, and often exhibits a decay in {v^ r) that resembles that of the variance of analytic white noise. 
If g\ {v^ t] is exactly constant across the a mb iguit y plane we can propose an estimator of a\ {v, r] given by equating o'xi^^^] 



to the estimator from \26\ . Propositions |3 .2| and |3 .3 argue that in regions of the ambiguity plane o'x{i'^ r] is smooth, and can be 
approximated as constant over a given region, again using (26). Noting that as the process is underspread t] | <C cr^(z/, r] 

for most (i^, r], and so thresholding is an admissible estimation procedure. We therefore propose a Local TEAF (LTEAF), by 
estimating the variance in given regions, just like the LBTEAF, but without removing the bias. We will divide the ambiguity 
plane into regions according to Fig. [TJa) to obtain both the LTEAF and LBTEAF. Using separate regions in the local time- 
frequency domain to estimate the variance of the EMAF is a natural procedure also for stochastic processes. It is reasonable to 
assume that the time-varying spectrum is smooth in global frequency and global time. If the time-varying spectrum is a member 
of a Sobolev space of order k, then the modulus of the Fourier transform must decay faster than zy~(^+^) in local frequency, 
and r~^^~^^'^ in local time, see for example |25|. It would therefore be natural to use regions defined in terms of distance from 
(ly^r) = (0,0) as Mk, but to more naturally meld with digital sampling we propose to use the set of square annuli, similar 
to sampling discussed in 1261 . Johnstone and Silverman in a similar spirit |[T2ll proposed level dependent thresholding using 
a different variance for each wavelet level, in 1-D. As a final note the proposed estimators may not correspond to valid AFs. 
Only AFs that are the FTs of valid autocovariance sequences are valid. We do not think this corresponds to a major flaw of the 
procedure, as we are mainly concerned with functions whose support will inform us of the resolution of the autocovariance. 



A. Estimated Spread 

Using thresholding methods, we have produced an estimate of the AF in the entire ambiguity plane. We also propose 
an estimator of the spread of the estimated AF, see O, ifTTl . Note that for processes that are not strictly underspread, our 
thresholding procedure will identify regions where the mean of the EMAF is sufficiently distinct in magnitude from the variance 
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of the EMAF. This in essence corresponds to comparing the magnitude of the AF at a point, with the magnitude of the AF 



at other points, see (32). We additionally note from (32) that if the AF is mainly centered at other {y' ^r'] then the variance 
of the EMAF will be large compared to the mean of the EMAF. In this instance the contribution at {v^ r] is not significantly 
contributing to the structure of X[t]. We define the ambiguity indicator cell to be 



a(^,r]=/(|Axx*(^,r]|^0), ve 



1 1 

2' 2 



r G Z. 



The total spread of the AF of X[t] over a time-frequency region S is given by 



< (5) = 



< 1. 



A process X[t] is strictly AF Compressible if {S) << 1. We define 



Cx{v,T]=li\Afl,{v,T]\^Q), 



1 1 

"2' 2 



,TG [-(7V-1),(7V-1)], 



and thus 



ix {S) 



(33) 



(34) 



(35) 



(36) 



is an estimator of the total spread of the process. Matz and Hlawatsch define extended underspread processes, as those whose 
spread is essentially limited. Such processes will be approximated by the TEAF to the region of their essential support, i.e., 
where their magnitude is non-negligible in comparison to the rest of the ambiguity plane. This permits the quantification of the 
degree of variability of the time series at each lag, and the band of variability supported at a given lag tq, can be determined 
by 6c : r = To}). 



V. Examples 

We estimate the mean squared error (MSE) by comparing the estimated AFs to a theoretically based quantity that is 
the expected value of the EMAF, for a finite value of N, limited to the support p of the AF of the process of interest, 
t]/ ((i/, r] G p), and refer to this as the A^-AF. If the process is the analytic process constructed from a real- valued 
process, then p is the support of the AF of the real-valued process. Note that Axx* {^^ t] is based on the infinite sequence 
{Mxx* [tj'^]}trez ^^^^ exhibit any finite sample issues with spreading in local frequency and time lag due to finite 
sample effects. We can never observe such values in a finite digital sample, because the maximum concentration of energy 
will behave like N (our number of sample points), and thus will be finite for finite N. Therefore, we ideally insist on the 
concentration of the AF to where Axx*{^^'^] is supported, but only letting the A^-AF take finite sample length realizable 
values. 



A. Linear chirp immersed in white noise. 

We first consider the case of a deterministic linear chirp, g[t\ = exp[j7r(2at + immersed in analytic white noise. 

Chirps are commonly characterised in the ambiguity plane |7|. The chirp has a starting frequency a = 0.1, with chirp rate 
P = 9.0196 • 10~^, and the noise has variance = 0.6. We generate data sets of length N = 256 from this process. We find 
the 7V-AF of the chirp 

^xi* ^] = exp {jTT [2aT - (3t^ + {(3t - v){N ^ |r | - 1)] } DN-\r\ {Pr - v) (37) 

for Pr = v and zero otherwise. We have not inserted jSr = v m this equation, since we are dealing with discrete values and 
Pt will not be identically equal to v anywhere. Instead, we will for each find ui = min \f3Ti — y\. Thus, the A/'-AF will be 
zero except for at the points {ui^Ti), where it is defined by (37 ). Fig. [2] (a)-(c) shows the EMAF, TEAF and LBTEAF for one 



realisation of the process. The EMAF is substantially corrupted by noise, but both thresholding methods have removed most 
of the noise. It should be noted that the line has exhibited some increasing thickness especially towards high |r|. This will 
explain why the MSE results are not quite as good as might be anticipated from this figure. To quantify the performance of the 
estimators, we will estimate the MSE by Monte Carlo simulation, where we compare with the A^-AF. We generate K = 5^ 000 
realisations of the process. The estimated MSE is shown in Fig. [2jd)-(f). Thresholding has reduced the MSE a great deal. 
To observe differences between the two thresholding procedures, note that the TEAF suffers from estimating a single noise 
variance, in that the regions where h{i'^r] inflates the variance, noise remains in the estimation (see for example the band 
around the chirp in Fig. [2je)). The LBTEAF looks like a considerable improvement, as it removes more noise, but does suffer 
from difficulties in estimating the noise around the rim of the domain. This leaves a "badly cleaned window" effect. 

To quantify our visual impression of these procedures, we look at the total MSE of each estimator, which is obtained 
by summing the MSE over the i^-r-plane, averaged over the five thousand realisations. The resulting total MSEs and their 
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(b) 
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(c) 
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(e) 



(f) 



Fig. 2. Chirp immersed in white noise, (a) EMAF (b) TEAF (c) LBTEAF for one realisation. Estimation of the MSE for (d) EMAF (e) 
TEAF (f) LBTEAF. All are on a dB scale. 

standard deviations are shown in Table |l] We see that the local bias corrected thresholding yields the lowest MSE, but both 
thresholding methods have quite significantly reduced the MSE of the EMAF. The reduction, is as mentioned, not quite as large 
as might have been anticipated, as the chirp has been broadened near its true support. But the reduction in MSE is respectable 
corresponding to a factor of almost four. 

We estimate the total spread for the chirp immersed in analytic white noise by Monte Carlo simulation. The estimated total 
spread and the standard deviation of the total spread is given in Table [ll| The EMAF will not be zero anywhere, so we do not 
need to estimate the total spread for this estimate, but equate this to 1. Theoretically, the AF of the chirp should be nonzero 
for only 511 of512x511 cells, which gives us a spread of 0.002. Our estimated spread is larger than this number, but reflects 
the sparsity of the TEAF. 



B. Stationary process 

We consider the analytic process X[t] corresponding to the real- valued stationary Moving Average (MA) process R[t] = 
X^ilo — where ^ ~ A/'(0, cr|) and E — t]} = cf'^5[t]. The autocorrelation of this process is 



L-\t\ 



for Irl < L. 



(38) 



i=0 



We generate reaUsations of the process with w = [l 0.33 0.266 0.2 0.133 0.066]^, length N = 256, L = 5, and 
cr| = 1. Fig. [3|a)-(c) show the EMAF of one realisation of the process and the result of thresholding the EMAF. Both 
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Chirp 


MA 


UM 


TVMA 


Total MSE of EMAF 


7.5382 • 10' 


2.0337 • 10^ 


3.3225 • 10' 


5.0842 • 10' 


Total MSE of TEAF 


1.9622 • 10^ 


9.7961 • 10^ 


1.7409 • 10^ 


3.1209 • 10^ 


Total MSE of LTEAF 




9.5598 • 10^ 


1.7076 • 10^ 


2.951 • 10^ 


Total MSE of LBTEAF 


1.8398 • 10 








STD of total MSE of EMAF 


9.9790 • 10^ 


4.6983 • 10' 


7.3221 • 10^ 


1.3809 • 10' 


STD of total MSE of TEAF 


4.0799 • 10^ 


4.9168 • 10^ 


5.1692-10^ 


1.7156-10^ 


STD of total MSE of LTEAF 




4.1841 • 10^ 


4.6430 • lO'* 


1.1522 • 10"^ 


STD of total MSE of LBTEAF 


4.0407 • lO'' 









TABLE I 

Average total MSE and the standard deviation (STD) of total MSE. 





Chirp 


MA 


UM 


TVMA 


Total spread of TEAF 


0.013 


0.0013 


6.5513-10-^ 


0.0012 


Total spread of LTEAF 




0.001 


6.3735 • 10"^ 


0.001 


Total spread of LBTEAF 


0.0156 








STD of total spread of TEAF 


0.0028 


0.0016 


8.9984 • 10-^ 


0.0018 


STD of total spread of LTEAF 




0.0012 


8.1374-10"^ 


0.0014 


STD of total spread of LBTEAF 


0.0039 









TABLE n 

AVERAGE TOTAL SPREAD AND STD OF TOTAL SPREAD. 



thresholding procedures have zeroed out most of the points for 7^ 0, and correspond to substantial improvements to the 
EMAF. Fig.|4|;a) shows the EMAF and the thresholded EMAF for one realisation with the 7V-AF for = and r G [-10, 10]. 
We see that the thresholded EMAF and the N-AF are very similar, and only when the variance of the EMAF become comparable 
to the modulus square of the AF is the process thresholded. This exactly corresponds to the support of the real-valued process' 
autocorrelation. 

The results of the MSE estimation for these estimators are shown in Fig. [3jd)-(f). Again, in most regions of the ambiguity 
plane we have a satisfying small value of the MSE, but at the ends, due to the low value of the threshold, suffer from the 
previously noted effects. Thresholding has reduced the MSE, and we see from Table |l| that the MSE is actually reduced by 
more than a factor of 100, which is a substantial reduction. Note also that the LTEAF has lower MSE than the TEAF. We 
estimate the total spread of this process, see Table [ll| The theoretical spread is in this case equal to 11 cells out of 512 x 511 
which is 4.2044 x 10~^. Our estimated spread is larger, due to the analytic signal spreading energy in r and the points at 
the ends that have not been successfully thresholded, but still reflects the geometry of the support. Using the discrete analytic 



signal and problems with edge effects makes the inflated spread inevitable. Whilst propositions 3.1-3.3 give the decay of the 
variance of the EMAF, at the rim of the ambiguity plane such expressions are not accurate as the 0(1) error term will be of 
comparable magnitude. 

C. Uniformly modulated process 

Next, we consider the analytic process X[t] corresponding to the real- valued, non- stationary. Uniformly Modulated (UM) 
process, R[t\ = ax[t]£,[t], where ^ ~ A/'(0, 1) and E{^[t]^[t — r]} = S[r]. The time- varying variance is defined as crj^[t] = 
sin^(27r/ot), which has an FT of S^x* (^) = I {2S{iy) - 5{y - 2/o) - 5{u + 2/o)) . The 7V-AF is A^xx^ {v, r] = for = 
andr = 0, A^xx*(v,t] = -7V(l/2 - |2/o|) for u = ±2/o,r = 0, and zero otherwise. We generate a realisation of the process 
of length = 256 with /o = 0.09. The EMAF and the results of the two thresholding procedures based on this realisation 
are shown in [5ja) - (c). The EMAF is very noisy, but both the TEAF and the LTEAF are substantially cleaner. In Fig. |4jb) 
we see the EMAF and TEAF for r = 0, and we observe that the thresholding has kept only the three points specified by 
the N-AF. Estimates of the MSE of these methods are shown in Fig. |5jd)-(f). Again we see that the thresholding has greatly 
reduced the MSE. From Table |T| we see that also for this process the MSE is reduced with a factor over one hundred from the 
EMAF to the TEAF and LTEAF, with the LTEAF doing a bit better than the TEAF. The theoretical spread of this process is 
1.1466 X 10~^, and our estimated spread is given in Table |ll| 

D. Time-varying MA 

Finally, we will combine the MA with the uniformly modulated process, thus definin g a re al- valued Time- Varying MA 



(TVMA) as R[t] = crx[t] Y^f=oU)i^[t — i], where the w's and ^[t] is as defined in Section V-B We use crx[t] = sin(27r/ot) 
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Fig. 3. MA process, (a) EMAF (b) TEAF (c) LTEAF for one realisation. Estimation of the MSE for (d) EMAF (e) TEAF (f) LTEAF. All 
are on a dB scale. 



with /o = 0.042, and we generate samples of length TV = 256. As always, we work with the analytic process corresponding 
to the real- valued process. We find the N-AF of this process as 

\r\) lo^' [S{f - fo) + S{f + fo)] ei^-S-df for v = QM<L 

^1/2-2/0 ^ /o)eJ-2-/-d/ for V = 2/o, |t| < L (39) 

for V = -2/0, |t| < L 



{N 

-(^-M)/o 



{N-\r\)ff'S{f-fo)e^^-f^df 



where S{f) is the FT of the autocorrelation in (38). We show the EMAF, TEAF and LTEAF of one realisation in Fig. |6[a)-(c), 
and we see that the thresholding has worked well, retaining only a few points. Fig. |4jc) shows the N-AF, the EMAF and the 
LTEAF for = 0. Likewise, Fig.[4jd) shows the EMAF and LTEAF at r = 0. These two figures demonstrate that the geometry 
of the non-stationary process is retained by the thresholding estimator, even if the spread of the EMAF is cut short when the 
variance of the EMAF becomes too large. The estimated MSEs are shown in [6|d) - (f), and the total MSE in Table [I| The 
MSE for the LTEAF shows a distinct improvement, especially near the region (0,0). The thresholding methods have again 
reduced the MSE with a factor over one hundred, and the total estimated spread demonstrates that the AF of this process is 
extremely sparse. 



VL Conclusions 

In this paper we have introduced a new class of estimators for the AF of a non- stationary process that exhibits sparsity in 
the ambiguity plane. The AF is a fundamental characterisation of a non- stationary process, and many important properties of a 
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Fig. 4. (a) TEAF (dashed), EMAF (dotted) and A^-AF (solid) for zy = for one realisation of the MA process, (b) TEAF (soUd) and 
EMAF (dotted) for one realisation of the UM process for r — 0. (c) LTEAF (dashed), EMAF (dotted) and N-AF (solid) for ly — for one 
realisation of the TVMA process, (d) LTEAF (solid) and EMAF (dotted) for one realisation of the TVMA process for r = 0. 



process can be deduced from its AF. The inherent resolution of the process in global time and frequency is an example of such. 
The AF has been used in estimating the generating mechanism of a non- stationary process |8|. It is also a popular tool in radar 
and sonar |5|. Despite this fact little efforts have been focused on the estimation of the AF, especially to produce estimators 
amenable to the determination of spread (size of the support). Characterisation of limited support is vital in designing the best 
analysis tools for generic second order non- stationary processes, a problem that remains open 1 1|. 

Based on the assumption of compression of the AF (small support), we proposed different threshold procedures for estimating 
the AF. The advantage of using the threshold estimator is that the size of the magnitude of the AF is compared with the magnitude 
of the AF at other time-frequency cells. Only if the local magnitude is sufficiently large is the value not thresholded. This 
should enforce a strict support from for example a process whose AF is only extended underspread, see |2|[p. 1072]. Unlike 
Matz and Hlawatsch we do not calculate the spread of the process by fitting the support of the AF into a box centered at 
(0,0), but simply count the number of non-zero AF coefficients. Conceptually this can be thought of as determining a finite 
number of cells that would represent most of the structure of the observed process. If a finite number of cells represent a 
full process then estimation of its generating mechanism is possible. Hence whilst it clearly is not necessary to constrain the 
process to be underspread, i.e. be concentrated in support around the origin in the AF domain |4|, it is necessary to constrain 
the possible degree of dependency in the process to ensure that the generating mechanism of the process could be estimated. 

Pivotal to thresholding procedures is determining a suitable threshold. We implemented the thresholding with a global 
variance estimate as well as a local variance estimate, both adjusted for each point in the ambiguity plane. For deterministic 
signals in analytic white noise we also proposed a procedure which removes bias in the estimator caused by the noise. We 
demonstrated the superior properties of the threshold estimator in finite sample sizes for a variety of common types of non- 
stationary and stationary processes. The reductions in MSE compared to a previously used, if naive, estimator for the AF 
were considerable, up to factors of over a hundred. Formally our proposed threshold procedure for zero-mean processes is 
only valid if the process is underspread in which case we determine the asymptotic distribution of the EMAF. We conjecture 
that asymptotic normality can be shown for a larger class of processes, and that in the case of any process with compressed 
AF, thresholding is an appropriate estimation procedure. The estimated spread does not fully reflect the degree of sparsity of 
the AF: an inflation is accrued due to the usage of the discrete analytic signal and edge effects. Resolving such effects fully 
remains an outstanding issue. The definition of spread in this paper has the clear advantage of interpretation, corresponding to 
the fraction of points where the mean of the EMAF dominates the variance of the EMAF. 

The AF remains the least studied of time-frequency representations of non- stationary signals, perhaps because its arguments 
lack direct global interpretability. A large class of processes that can be estimated exhibit compression in this domain. We 
anticipate that further study of the AF of non- stationary processes will yield understanding into what generating mechanisms 
can be determined from a given sample with fixed sampling rates and sample size. 
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(e) (f) 

Fig. 5. UM process, (a) EMAF (b) TEAF (c) LTEAF for one realisation. Estimation of the MSE for (d) EMAF (e) TEAF (f) LTEAF. All 
are on a dB scale. 



Appendix A 
Proof of Proposition s. II 

The Empirical Cross Ambiguity Function for two finite equal length samples of signals and {^2[^]}t of length 

is given by A,,,>^{v,t] = Y.^IITS)'"'' aMd^it " rje-^'^-*. Then 

Ar-l+min(0,r) Ar-l+min(0,r) 

Axx'{v,t] = Yl Mxx'[t.r]e-i^^-* = ^ ^g[t] + W[t]){g*[t- T] + W*[t- t]) + e-^^^^* 

^=max(0,r) t=max(0,r) 

= Agg^ {v, t] + Awg- {i^, t] + Agw- {i^ , t] + Aw w * r]. (A-1) 



DEPARTMENT OF STATISTICAL SCIENCE RESEARCH REPORT 293 



14 




As g[t] is deterministic, and W[t\ is zero-mean, it follows E j^vt^^* (i^, r] | = 0, and E | A^vi^* (i^, r] | = 0. Using (4) for 
r > and the fact that Sxx* (^^ /) = ^vf^(^) for / > and zero otherwise for the noise process, we find for r > 



^0 

1/2 nl/2 



JO 

j-K: 



gj^.'(iV-l-r)^2^^(^/)^j2^(.'+/)r^^_^(^/ _ j^^^-jnu{N+r-l) 
.1/2 



^0 



df 



[w/g-J-KA^+r-l)2)^_^(^)gjW2ginc(^r/2), 



where (as usual) sinc(x) = sin(7rx)/7rx. Thus taking expectations of \A-\\ and using the linearity of E{-}, the result follows. 
Mutatis mutandis the c alculations are implemented for r < 0. 

We start from ( A-1 ) and note that var | Axx* (^^ | can be written in terms of covariances of Agg* (^^ AgW*{i^^'^] etc. 
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Here, g[t] is deterministic and W[t] is Gaussian proper and so 

var |A^^*(i/,r]| = 0, cov | A^^* (i^, r], A^^^* (i/, r] | = 

COY ^^Agg*{u,r],Agw*{^,r]^ = 0, coy ^^Agg{u,r], Aw w*{^,r]^ =0 

COY i^Awg*{j^,r],Agw*{j^,r]^ = 0, coy Aw g*{j^,r], Aw w*{j^,r]^ =0 

COY i^Agw*{j^,r],Aww*{j^,r]^ = 0. (A-2) 

Ar-l-|r| 

We next define G(/, r] = 9[t — I'^IH'^ < 0)]e~•^^^•^^ where /(•) is the indicator function, and note that G'(/, r] is 

t=o 

supported on negative frequencies as well even if g[t] is analytic. If N — \r\ oo, G(/, r] will tend toward only being 
supported on positive frequencies. Then, for r > 0, 

'N-l-T 

r) 



var{lH-^*(^,r]} = var i ^ W[t ^ r]g''[t]e-^^^^^'+^^ 

( ^ /N-l-r \* 



j27rr(/i-/2) 



= var <^ dX(/i)G*(/i - ^, rje^'^^^^^ #2 

= 1'^' 1'^' ^ {dX(/i)dX*(A)} G*(/i - z., r]G(/2 - z., r]e 
= / \G{f-iy,r]\^df (A-3) 

After a change of variables f = f — are given an outer integral of J^^^ ^ over rather than Jq^^^. For z^ > we can 
rewrite this as J^^'^ ^ plus a contribution that will have the inner integrals integrating to negligible contributions. For u < 
we can rewrite this as J^^^ plus a contribution that will have the inner integrals integrating to negligible contributions. We 
denote h{u^r] = /mlx(^^o)^^^ \G{f^^]f • ^ generic harmonizable process r > gives 

_ . /-l/S rl/2 .1/2 nl/2 

var{Axx*(^,r]}= e.2-(/i-«i)rgM/i-/2-«i+a.)(iV-i-r)^^_^(^^ _ _ ^) 

This follows by direct calculation starting from (|4]). Using Isserlis' theorem 1271 we write the variance as 

nl/2 .1/2 .1/2 .1/2 

var j = / j j i Sxx-{h - oli,oli)Sxx-{ol2 - f2, f2) 

gJ7r(/i-/2-«i+«2)(iV-l-r)gi27r(/i-ai)r^^_^(y.^ _ _ j,)Dn-Mi - 0L2 - v)dhdf2dOLxdOL2. (A-4) 

which for analytic white noise reduces to 

var |Axx*(z^,r] j = \ J ^^-r(/i - /2 - ^)^/i#2 

r>l/2+min(0,-zy) ng2{N-T) 



<^w{N-t) / sinc2(e)d^d52 

Jmax(0,-z^) J-{l/2-g2){N-T) 

<T^(iV-r)(l/2-|H) + 0(l). 



Implementing the same calculations for r < derives analogous expressions with r replaced by |r|. Putting (A-2), (A-3) and 
var {Agw* = var {A^^^^* (— z^, — r]} together, and mutatis mutandis implementing the calculations for r < 0, thus yields 

var{lxx*(^,r]} = + a^/i(z., r] + a^/i(-z., -r] + a^(Ar - |r|)(l/2 - |z/|) + O (1) . 
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We Start from (A-1) and note that rel 't] | can be written as an aggregation of relations. As g[t] is deterministic 

and W[t] is zero-mean Gaussian, 

rel{l^^*(i/,r]} = 0, rel {l^^* (i/, r], A^y^* (i/, r] } =0 (A-5) 
Te\i^Awg*{j^,r],Aww*{j^^r]^ 0, rel | A^h^* (i^, r], A^^h^* (i/, r] | =0. 



For r > 



rel ^Awg* {y, t] , Agw* {y, ^] } = E \Awg- (z^, r]Agw* (z^, r] | 

r .1/2 .1/2 _ _ ^ 

.1/2 

= / - v,T]G{f + ^, -r]e^-^-^^-^-)-4f. (A-6) 

The factor e"^"^^^^ is not present for r < 0. We define h'{v,T] = G*{f,r]G{f + -T]e^'^''^f-''^<^>^^ df , 

where sign(r) = 1 for r > and sign(r) = — 1 for r < 0. If G(/, r] is mainly supported only over a range of frequencies we 
would assume with increasing v that h'{v^ r] 0. Furthermore we note that rel | A^y^* {v^ r] | =0, and rel | A^^f* (z^, | = 0- 
We can note that for a generic harmonizable process X[t\, 



.1/2 .1/2 .1/2 .1/2 

.■e.{.4.v.(.,.l}=/ III 

Dn-tUi -f2-y) DN-r{ai -a2- iy)E {dX(/i)dX*(/2)c/X(ai)dX*(a2)} - E {l^x* (z^, r] }' . 

This follows by direct calculation starting from (|4]). Using Isserlis' theorem and duplicating the calculations for the variance, 
we write the relation as 

.1/2 .1/2 .1/2 .1/2 

Te\\Axx*{j^,r]\= / / / (/i - Q^2, Q^2)5'xx* (o^i - /2, /2) 

^ Jo Jo Jo Jo 



gJ>(/i-/2+ai-a2)(Ar-l-r)gj27r(/i+ai)rg-2j7ri/(Ar+r-l) 



Jvr(/i-/2+ai-a2)(iV-l-r)gi27r(/i+ai)r^^_^(y.^ _ _ _ _ jy)e-^^''^^^^^-^^ dfidf2daida2. (A-7) 



For an analytic white-noise process this corresponds to 

rl/2 .1/2 .1/2 .1/2 



_ . ri/-^ ri-/-^ pi/-^ rW^ 

Te\{Axx*{iy.r]] = J J J J ^^^(/i - c^2)a^^K - /2)e^-(^-^^+— 
g,2.(/,+a,)r^^_^(j^ - /2 - i^)DN-r{ai - ^2 - u)e-^^^^^''^^-'Uhdf2daida2 

= f'^' /'^'a^e^--(^^-^^+^^-^^)(^-^--)e^-2-(^^+^^)-I)^_,(/i-/2-z.) 
Jo Jo 

DN-r{f2 - /i - ;^)e-2^'^-(^+--i)d/id/2. 



If u = O {1/{N — r)) then this corresponds to (up to order {N — r)/N) 



l- J Jo J(f-i/2)(N-T) \(N-t) J \ (N - t) J N-T 

70 A.f 



'^W ~ • 

Jo J(f-l/2-u)(N-T) N-T 



If ly is not O {1/{N - r)), then rel |lxx*(z^,r]| is negligible. If u = O {1/{N - r)) then with ^(z/) = (A/" - r)z/ = O (1) 
we obtain that 

r»l/2+min(0,zy) 



rel { Axx* (^^, r] } = (iV - r) / sine (/') sine (/' + 2^(j.))e-2i--(^+-i) d/' / 



pjf47r(l/2+min(0,i^))r pj47r max(0,z/)r 

= (iV - T)a^L {N - T, v) e-2j--(iv+r-i) ^ ^ f ^ ^ 



-(iV-T)<j^L(iV-T,i.)e 



^H^(|l;M,(,^0)-i/(r = 0) 



■0(1). 
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This defines the quantity L {N — r^u), which decays Hke O ^^^^^ Thus if r 7^ 

rel [Axx^ r] ] = -{N - r)\u\a^L (TV - r, u) e-^^^^^^-'^smc{2\iy\T) + 2a^^h\iy, r] + O (1) . 
Mutatis mutandis calculations can be replicated for r < 0. 

Appendix B 
Proof of Proposition [3T2l 

For r > we note that 

Note that DN-r{^) is an even function, and the expression is valid for v < 0. Calculations can mutatis mutandis be replicated 
for negative r. We note that with r > 

var{Axx*(^,r]} = J J ^xx* (/2)l^^-.(/i - /2 - ^M/i#2 

= (7V-r)(l/2-|^|)Axx*(-^,0]+O(l), 

with Axx*{j^,^] = Imax{o^-u) '^'' ~ ^) S X X* {f)^^'^^ df/ {1/2 — Mutatis mutandis the calculations can be 

replicated for r < 0. Finally we start from \A-1) and note that for stationary processes 



g jT:v{N+T-l) 



rel{Axx*(^,r]| = J J ^xx* (/2)e^'"^^^+^^^" 

DN-r{fl - /2 - ^)e-2^--^(^+--l)l)^_,(/l - /2 - ^)6i/l#2 

/.I/2 />(Ar-r)(fir2-z^) / ^ \ 

= e-2i-(^+-i) / / 5xx.(52)5xx. {92 - ^ - A e^^<'^^-^-^->^ 

Jo J{N-T){g2-l/2-v) \ J 

= e-2^-'^(^+--i)(iV - r)(l/2 - |r.|)i(iV - r, r] + O (1) . 

mutatis mutandis results may be derived for r < 0. 

Appendix C 
Proof of Proposition [33] 

Using ^ and by direct calculation we note that 

^1/2 nl/2 



-j7ru{N-\-T-l) 



Jo Jf-l/2 

/>l/2-max(z.,0) rfiN-r) / 1/ \ / 1/ \ 

Jmax(0,-z/) ^(/-l/2)(Ar-r) " ^ / V^^ " ^/ 

/>l/2 — max(z/,0) /•cxd 

^maxfO, — z/) J — 00 



I max(0, — z/) 

pj27r(l/2 — max(i/,0))r pjf27r max ( — 1^,0)) r 

Sxx.(^) + O (1) 

(1/2 - |i.|)Sxx.(^^)e^'"(i/2-.)rgj^^((i/2 - |j.|)r) + O (1) . 
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Mutatis mutandis we may derive the results for r < 0. We start from ( |A-4| ), and again by direct calculation 

.1/2 .1/2 .1/2 .1/2 

var{Axx*(^,r]}= / / / / ^xx* (/i - c^i, ai)^xx* (/2 - c.2, c.2)e^''"^^^-"^^^ 

^ Jo Jo Jo Jo 
gi-(/i-/2-ax+a.)(iv-i-x)^^_^(^^ ^)DN-r{ai - - v)dhdhda^da2 

fill 1-1/2 ff ra 

Jo Jo Jf-l/2Ja-l/2 

gj27r(/-a)r^^_^^^/ _ _ y)dv'da'dfda 

rl/2-v nl/2-v rf{N-T) paiN-r) / / _ / 

/ / / ^xx*{f-a)^^^^Jf--- a 

J-v J{f-l/2){N-T) J{a-1/2){N-T) \ 1^ - T 

^Mu'-a')iN-l-r)/iN-r)^j2nif-a)r-^^_^^^y^j^ _ r))D n -r{oi' / {N - T))dv' da' df da / {N - rf 

.l/2-max(iy,0) .l/2-max(zy,0) 
^max(0, — zv) ^max(0, — z^) 

= / [--\v\]\Y.xx^{f)\\^'-y-df^O{l) 
J-i/2+\y\ / 

= (1/2 - \v\){N - r) / I XX U I ^,2./.^^ ^ o ^ 

J-l/2+|i.| (^-^) 



We start from ( |A-7| ) 

r>l/2 .1/2 .1/2 .1/2 



.i/Z .i/Z .i/Z 

Te\lAxx*{j^,r]j= / / / / (/i - Q^2, Q^2)5'xx* (o^i - /2, /2)e' 

><^,2.(/,+a,)r^^_^(y^ - /2 - iy)DN-riai - a^ - v)e-^^^-^''^--^Uhdhda^da2 

.1/2 .1/2 .1/2 .1/2 

JO Jo Jo Jo 

eJ^-ih+-^)rDN_,{f, - /2 - iy)DN-r{ai - a^ - iy)e-^^^^^^^--'Uhdf2darda2 
We implement the change of variables u'' = {v' — v){N — r) and a" = {a' — v){N — r). 

.1/2 .1/2 n{f-v){N-T) n{a-v){N-T) / 

rel Axx*(^,r] = / / / / ^xxAf-a 

^ ^ Jo Jo J{f-l/2-v){N-T) J{a-l/2-v){N-T) \ 



j7^ifi-f2-\-ai-a2)iN-l-T) 



N-r 



7V- 



.l/2+min(0,z^) .l/2+min(0,i/) 

= e-^^^"- / / Sxx. (/ - a + (/ - - a)e^'2-(/+")- d/da + 0(1). 

^max(0,i^) ^max(0,i^) 

Mutatis mutandis the expressions may be derived for r < 0. 

Appendix D 
Proof of Theorem 13. II 

We outline the proof for r > 0, the same arguments hold mutatis mutandis for r < 0. We let be the analytic 

extension of a real- valued process or equivalently X[t] = Y[t] -\- jU[t]. As is a real- valued underspread 

process there exists a constant T = 0(1) such that V |r| > T, coy {Y[t],Y[t ^ r]} = 0. This implies that V |r| > T 
COY {Y[t]^U[t -\- r]} = 0(l/(r — T)), where is the discrete Hilbert transform of For convenience assume 

N — r = Tn. This is of course not a necessity but simplifies the proofs. Otherwise for Tm such that Tm < N — r and 
T(m + 1) > N — T we can split the sum in the EAF into two parts, and show the latter part has negligible contributions to 
the total sum. The first part, summing from zero over t up to Tm, can be treated as if A/" — r = Tm with m = n. We now 
take a full length observation Y = [Fq, • • • , Y^-i] and construct n sub vectors by the following mechanism 

Yi = [Y[0],Y[T], Y[{n - 1)T]] , Y2 = [Y[l],Y[T + 1], . . . , Y[{n - 1)T + 1]] , .... 
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The vectors are constructed so that the elements of are uncorrelated and if we define 11^ and X^, as the obvious 
extensions to Y^, then coy {Yk[t],Ui[s]) = O {l/\k - I ^ T{t - s)\), cov (X^M, X^Js]) = O {l/\k - I ^ Tit - s)\) if k ^ 1. 
Before constructing the CLT type of argument let us study how we shall represent Axx*{i^-,t]. We write 



T-l n-1 



The expectation of t] follows directly from this expression, and we additionally note that e|Mxx*[^,t]| 

Mxx*[t,T]. 

Using Isserlis' formula 1271 and the assumption that the process is proper, we find the variance of Axx* {^^ t] as 



-1 n-l T-l T-l 



Vi=0 V2=0 Ui=0 U2=0 

where 

CN{iy, r; ui,U2,vi,V2] = e-^'^^'^^^^^-^^^+^^-^^^M^x* ^iT + ^xi + r, T{vi - V2) + iii - U2] 

Mxx* i^iT + uuT{vi - V2) -^ui- U2]. 

As is exactly underspread if Vi / V2 and / ^2 i 1 we have that with ti = viT + ui, CAr(z^, r; ixi, 1^2, ^i, ^2] 

O (l/(ti -12^) . Furthermore, E|^i-^2|>i ^' ^i' ^2, ^1, ^2] = O (log(n)) . Otherwise 

cn{-i^,t]Ui,U2,vi,V2] = M^x* [^1^ + ^1 ^r,ui-U2^ {vi - V2)T]M^x*bi^ + ^i^^i - ^2 + (^1 - ^2)^]. 
Thus we find that (up to O (log(n))) 

T-l ni + (T-l) n-l 

var [Axx* ^] } = E E E e'^^^"^'''''''^ Mxx^ [vT ^ui^t.ui- U2]M*xx* [vT ^u^m- U2]. 

ni=0 U2=ui-{T-1) v=0 

Let X = vT -\- ui and r' = ui — U2. Then this expression is rewritten for r > as 

N-T-l T-l 

var{lxx*(^,r]} = ^ E e-^''"'^ M^x* [x + r, r^]Mi^* [x, r'] + O (log(n)) . 

x=0 r' = -(T-l) 

Thus var | A^x* r] | = O (X). We note that 

n-l n-l T-l T-l 

rel|Axx*(i^,r]| = X] X] X] ^ rAr(i^, r; ixi, 2x2, ^1, ^2](i^, r], 

vi=0 ^2=0 ni=0 n2=0 

rN{i^, r; wi, ua, i-i, 1-2] = e-^'2^'^(^("i+"^)+"i+"^+2")Mxx. [uiT + ui + t, T(ui - V2) + wi - W2 + r] 

X M^;^ . [wiT + wi , T(wi - W2) + ui - U2 - r] 

If \ti -t2+T\>T and |ii - - t| > T 

r.(.,r; = O ( _ + , '(,^ _ _ 

whilst if \ti - ^2 + r| > T and \ti - t2 - r\ < T 

rN{j^,r;ui,U2,vi,V2] =0 



ti-t2^r 

Thus if |r| > T, J]^^^^^ rAr(i/, r; ixi, ^^2, ^i, ^2] = O (log(n)) . If |r| < T E|^i-^2|>i ^' '^i' ^2, ^1, ^2] = O (log(n)) . 

If r < T and |vi - ^2] < 1 we have that rAr(i^, r; i/i, 1x2, ^i, ^2] = e-^'2^^(^(^i+^2)+ni+n2+2r)^^^^ ^^^j^ + ^xi + r, T(vi - 
'^2) + '^1 — '^2 + ^]^xx* + '^1, ^('^1 — '^2) + '^1 — '^2 — 't]. We deduce that with x = vT -\- ui and r' = Ui — U2 that the 
relation of Axx*{i^-,t] is 

_ r 0(log(n)) if |r| >T 
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This completes the proof of the order structure of the EAF of an underspread process. 

To prove a CLT for the EAF we add the extra assumptions stated in the theorem. We need to use results for random variables 
that are both non- stationary and correlated to determine large same results. We note that for r > 

N-T-l N-T-1 

Axx'{'^,t]= Y1 QN-\r\[t]+ J2 e-i'^''^*+^'>Mxx'[t + r,T]. 

t=0 t=0 

We note that {Q7v-|r| [t]^ t = 0, . . . , — r — l} is a triangular array of centred random variables, that these by assumption 
are strongly mixing, and have finite second moments. The constraint Peligrad |28 | makes on the correlation is satisfied because 
X[t] is strictly underspread, and the Hilbert transform only induces suitably decaying correlation from such a process. We 
have assumed the Lindeberg condition holds, and note ([17]) from our assumptions. Hence from Theorem 2.1 of Peligrad 1281 
the theorem follows. The condition in ([TT]) does not have to be stated for the real and imaginary part separately as the relation 
divided by the variance goes to zero for (i^, r) ^ V. We can deduce the exact form of the mean, variance and correlation from 
the previous part of the theorem. 
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